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ABSTRACT 

Kinetic approaches provide an effective description of the process of particle 
acceleration at shock fronts and allow to take into account the dynamical re- 
action of the accelerated particles as well as the amplification of the turbulent 
magnetic field as due to streaming instability. The latter does in turn affect the 
maximum achievable momentum and thereby the acceleration process itself, 
in a chain of causality which is typical of non-linear systems. Different kinetic 
approaches are characterized by different levels and types of approximations 
that also imply different computational times. Here we present the results 
of two such approaches: one which is mathematically rigorous but rather de- 
manding from the point of view of computational time, and the other which 
is computationally very fast but based on an ansatz that, while physically 
justified, is not rigorous. The identification of possible differences can be cru- 
cial in assessing the possibility of implementation of one such calculation in 
hydrodynamical codes for supernova explosions. Special emphasis is given to 
a discussion of the appearance of multiple solutions in both approaches. 
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1 INTRODUCTION 



Recent observations of non-thermal radiation from shell-type supernova remnants (SNRs) 
are showing that effective particle acceleration takes place in these astrophysical objects (see 



FunkI (120071 ) for a recent review). In particular the detection of spatially resolved non-thermal 
X-rays from thin regions close to the forward shock in several SNRs has showed the first clear 
evidence for strong magnetic field amplification as could be expected if the shock is effi c iently 
acceler at ing cosmic r ays. The recent detection of TeV gamma rays ( lAharonian et al.l ( l2004l . 
2006j |bl. I2OO5 , 2006c )) is likely to provide further information on the magnetic field amplifi- 
cation by adding to multifrequency observations and by possibly allowing the measurement 
of the shape of the TeV gamma ray spectrum. 

In the scenario with large magnetic field, low fluxes of gamma radiation through inverse 
Compton scattering (ICS) of electrons are expected, thereby leaving a hadronic origin of the 
gamma ray emission as a more likely possibility. This latter conclusion is however of more 
limited strength and requires further confirmation through a continuous effort to detect 
SNRs in gamma rays, not only in the TeV region but also in the GeV energy range, as 
should become possible after the upcoming launch of the GLAST gamma ray telescope. 

X-rays observed in the 1-10 keV energy range from SNRs are generated by synchrotron 
emission of relativistic electrons and the thickness of the brightness profiles allows one to 



estimate the total field in the s 



lock vicinity to be of order 100 - 500 u G in basically all cases 



in which measurements exist ( jVolk. Berezhko fc Ksenofontovl ( 120051 )). This appears to be 



the strongest evidence for efficient cosmic ray acceleration at SNR shocks, where here we 
use the word efficient to indicate that an appreciable fraction of the kinetic pressure at the 
shock is transformed into accelerated particles, which in turn change the dynamics of the 
shock. In this r egime the standard test particle theory (TPT) fails and a non-linear theory 
is required (see iMalkov fc Druryl ( I2OOII ) for a review). From the phenomenological point of 
view the introduction of the non-linear effects is crucial. 

The reaction of the accelerated particles onto the shock structure has bee n calculated 
within different approaches, both semi-analytical (see iMalkov fc Drurvi (20011) for a recent 
review) and numerical (e.g. Ellison. Baring fc Jones ( IQQgI ): Kang fc JoneJ ( 2006 ) and ref- 
erences therein). Semi-analytical kinetic approaches provide a detailed description of the 
process of particle acceleration at cosmic ray modified shocks, including particle spectra. 
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modified shock dynamics and ther modynamics, and have rece ntly been generahzed to in- 
clude magnetic field amplification ( lAmato fc Blasil (120051 . l2006l ) ) . 

The main limitation of these approaches is the fact that they all solve the stationary 
problem of acceleration. This problem is simply ill defined in the case of acceleration of 
protons (or nuclei) in sources such as supernova remnants: in fact in these cases there is 
no appreciable energy loss process, and stationarity cannot be reached. This is in principle 
true even in test particle approaches, but there the problem is limited to momenta close to 
the maximum momentum, while the shape of the spectrum at lower momenta is rather well 
defined and stationary. In the non-linear regime, the entire particle spectrum changes as a 
consequence of the time dependence of the maximum momentum, therefore the assumption 
of stationarity is in principle not well justified, though it is tacitly assumed that the situation 
to describe is one of quasi- stationarity. 

Several approaches t o time-dep e ndent shock acceleration a t cosmic ray mo d ified shocks 



have been put forward (IBelll (119871 ): iFalle and Giddingsl (119871 ): iKang fc .Tonesl (120061 )). but 
all of them are based on numerical a pproa ches (e.g. finite difference solution of the equations). 
In particular, iFalle and Giddingsl (119871 ) found that quasi-stationary solutions are rather 
good approximations to the full time-dependent solutions. Despite this, a generalization 
of the semi-analytical approaches to include time dependence would be highly desirable, 
and would probably help to shed some light on the problem of the appearance of multiple 
solutions, as discussed below. 

The first kinetic mo del that p rovide d a s emi-analytical (stationary) s olutio n of the prob- 
lem was developed by iMalkovl (119971 ) and iMalkov. Diamond fc Volkl (120001 ) for strongly 
modified shocks with a given, spatially const ant, d if fusion coeffi cient. A simple, approxi- 
mate kinetic approach was later proposed by BlasJ (2002, 2004). This approach provides 
an accurate description of the shock modification and its effects on particle acceleration 
provided the diffusion coefficient has a sufficiently strong dependence on momentun i p (the 



Eichler 



under lying assumption is somewhat similar to the starting assumption adopted by 

mm)- 

It also allows to obtain the solution to the problem in a very short computational time, 
which is particularly important when the particle acceleration process must be described 
in the context of com plex and time consuming hydro dynamical simulations (for instance 
in the case of SNRs). iBlasi. Gabici fc Vannonil (120051 ) showed, among other things, that 
the solutions obtained with Malkov's method and those of Blasi's method are in very good 
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agreement for Bohm diffusion coefficient, while the agreement becomes gradually worse 
for Kraichnan {D{E) oc E^^"^) and Kolmogorov {D(E) oc E^^^) diffusion coefficients, as 
expected. Even in th e se tw o cases, however, the differences were not very large and the 
model of iBlasil (120021 . 120041 ) still provided an acceptable description of the main physical 



aspects of the problem. 

More recently a solution of the system of equations des cribing the transport o f accelerated 
particles and the dynamics of the shock was found by lAmato fc Blasil (120051 ). where the 
diffusion coefficient could formally be taken as an arbitrary function of both momentum 
and spatial coordinate (although solutions were actual l y corn puted for a few specific choices 
of D{x,p)). Contrary to the simple approach of BlasJ (2002) where a physically reasonable 
ansatz was used to simplify the problem, this approach is mathematically rigorous but more 
expensive from the compu tational point o f view. 

In a second paper by lAmato fc Blasil (120061 ) the self-generation of the waves and the 
determination of the resulting diffusion coefficient, though in the context of quasi-linear 
theory, were introduced. The method for taking into account magnetic field amplifica- 
tion during the acceleration process at a modified shock, as assessed in this work, also 
allowed for the determina t ion of the maximum momentum in these complex circumstances 
( Blasi. Amato fc CaprioU (20o3)). 

Since the two methods introduced above aim at describing the same problem but using 
a different set of approximations and, very important, with a quite different computational 
cost, we wish to assess here the advantages and disadvantages of using one or the other. 
While doing so, we also discuss an important feature which is found in both approaches, 
and mo re generally in semi-analyt ical kinetic approaches (as well as in stationary two-fluid 



models 



Drury fc Void (jl980l . Il98ll )). namely the appearance of multiple solutions. 



We conclude that for most applications related to the descr i ption of the phenomenology of 



particle acceleration in supernova remnants the model of [Blasil (|2002| . |2004J ) can be succesfully 
applied despite its simplifications, provided the diffusion coefficient scales rapidly enough 
with momentum. This is certainly the case for Bohm diffusion. Moreover we complete this 
simple model with recipes that allow us to determine the level of magnetic field amplification 
and the maximum energy of the accelerated particles. 

The presen t pap e r is s tructured as follows: in ^ we briefly introduce t he two kinetic 
approaches of iBlasil (120021 ). model A hereafter, and lAmato &: Blasil (120051 ). model B. In 
the subsections we discuss in detail the results of the two theoretical approaches. In ^ we 
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describe the introduction in the first approach of a recipe for the acceleration time which 
allows us to estimate the maximum momentum of the accelerated particles. We conclude in 

m 



2 TWO KINETIC APPROACHES 



In this section we briefly discuss two kinetic approaches, the one of lBlasil (120021 . 12004 ) (Model 
A) that is expected to work for diffusion coefficients D(p) that change rap idly enough with 
the particle momentum, and the one proposed by lAmato fc Blasil (120051 ) (Model B) and 
valid for any choice of D{p,x). The main point of the section is to illustrate the results 
obtained with the two approaches, identifying the advantages and possible shortcomings of 
the first, which is being implemented in some computationally heavy hydrodynamical codes, 
in which the formal solution of Model B would be hardly possible to introduce due to the 
much longer computational time required. 

The basic equations are the transport equation: 



d_ 

dx 



d 

D{x,p)—f{x,p) 



^ df{x,p) ^ 1 

dx 3 \dx , 



and the conservation equation for the total momentum: 
1 . 1 



Ux) = 1 + 



U{x) 



;U{X) 



-79 



(1) 



(2) 



'0 iM 

Here f{x,p) is the particle distribution function in the shock frame, u{x) is the velocity of 
the background fluid, which equals U2 downstream and changes continuously upstream, from 
Ml immediately upstream of the subshock to uq at upstream infinity. The quantity U{x) = 
u{x)/uq is the normalized velocity, bound to equal unity at a; — oo. In Eq. |2l C,c{x) = 
Pcr{x)/ pquI is the cosmic ray pressure at the position x, normalized to the kinetic pressure 
PquI at upstream infinity. It is customary to introduce the compression factor Rgub = Ui/u2 
at the subshock and the total compression factor Rfot = Uq/u2- Assuming homogenisation 
of the cosmic ray plasma in the downstream sectio n (df / dx~ = 0, a consequence of the 
assumption of stationarity) one easily obtains (IBlasil (120021 )) that the distribution function 
of the particles at the shock location is 



/o(p) 



3R, 



tot 



Wo 



RtotUpip)-! AttpI 



exp < - 



p dp' SRtotUpjp') 

Pinj P' RtotUp{p') — 1 



where we introduced the function Up{p) = Up/uo, with 



Mr 



Ml 



dx{du/dx)f{x,p) . 



(3) 



(4) 
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This result is very general and is not related to the assumptions of the model of lBlasil (120021 ) . 
Here we assumed, as usual, that the injection is a delta function in space (at the shock) and 
in momentum (at the injection momentum pj„j) in the form Q{x,p) = ^'^^""j^"^ 6{p—pinj)S{x), 

^ inj 

with Ugas,! = noRtot/ Rsub the gas density immediately upstream (x = 0~) and rj the fraction 
of the particles crossing the shock which are going to take part in the acceleration process. 

The compression factors Rsub and Rtot, if the gas in the upstream frame evolves adiabat- 
ically, are related through the following expression: 



R 



tot 



(7, + ml, - (7. - 1)^: 



79 + 1 

sub 



(5) 



Model A 



This model was proposed bylBlasil (l2002l ) and put in a form to include the reacceleration 
of seed particles by iBlasil (120041 ) . It is characterized by a simple implementation that makes 
it fast in terms of numerical computation while keeping all the main physical ingredients of 
the problem. For these reasons it has been implemented in hydrodynamical codes in order 
to cal culate the ins t antan eous spectrum of accelerated particles in an expanding supernova 
shell (lEUison et all (boOTI)). 

The essence of the approach is based on the observation that the function /(x, p) in Eq. H] 
is expected to suffer an exponential suppression at a distance Xp such that dxu{x) / D{p) ~ 
1. The spatial position where this suppression is found in general depends on momentum if 
the diffusion coefficient depends on momentum. These spatial locations for different values 
of p are w e ll def ined if D(j)) is a strong function of momentum, as was first pointed out 
by Eichler (IQtJ, and become more ill defined for Dip) weakly dependent on p. As a first 



approximation we may assume that the distribution function is f{x,p) = fo{p) at \x\ ^ \xp\ 
and vanishes at |a;| > \xp\. With this assumption Eq. Hlgives Up ~ u{xp). The distance Xp has 
the meaning of the typical distance to which particles with momentum p can diffuse away 



from the shock in the upstream fluid, and the fluid velocity at that point is u{xp) ~ Up. It 
follows that Eq. [2] can be transformed into an equation which depends only on the particle 
momentum p: 

1 1 



Up) = 1 + 



2 P 



(6) 



where 
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47r /"P" 



'max 



Up) ~ o — 2 / dpp v{p)fo{p), (7) 



and v{p) is the velocity of particles with momentum p. Differentiating Eq. [6] with respect to 
momentum we finally have: 

477 



dU.p 
dp 



^/•/-(79+l) 



„ ^pMp)fo{p). (8) 



'■0 

For p pinj we have Up — > ui/uq = Rgub/Rtot, while for p — »• Pmax one has f/p ^ 1. The 
procedure to calculate the solution is therefore straightforward: for a given guess of Rsub (and 
therefore of Rtot through Eq. [5]) we can solve the differential equation Eq. [8] coupled with the 
expression for fo{p) as a function of Up (Eq. [3]) in an iterative way. In general we end up with 
Up{pmax) 7^ 1; whlch Impllcs that the given value of Rsub is not a solution of our problem. 
The solution corresponds to that value of Rsub (and Rtot) for which Up{pmax) = 1- This value 
also identifies completely the function Up and the distribution function at the shock fo{p), 
which is the sought after result. It is important to notice that this procedure does not lead 
only to the determination of fo{p) but also to the values of the thermodynamical quantities 
at the shock, such as the temperature of the gas upstream and downstream. What Model 
A cannot provide is the spatial dependence of the quantities in the precursor, although we 
provide below a physical ansatz that may indirectly make this piece of information available. 
We later check the results by comparing them with the outcome of Model B. 

In b oth approaches that we de s cribe h ere we adopt the inje ction recipe known as thermal 
leakage f Blasi. Gabici fc Vannoni ( 2005 ): Gieseler et al.l (2000)) which connects the value of 
rj to that of Rsub in a unique way through the relation: 

V=^iRsub-l)ee-^'. (9) 

Here ^ is a parameter that identifies the injection momentum as a multiple of the momentum 
of the thermal particles in the downstream section {pmj = ^Pth,2)- The latter is an output of 
the non-linear calculation, since we solve exactly the modified Rankine-Hugoniot relations 
together with the cosmic rays' transport equation. For the numerical calculations that follow 
we typically use values of C, between 3 and 4, which correspond to fractions of order ~ 
10~^ — 10^^ of the particles crossing the shock to be injected in the accelerator. 



Model B 



A formal solution of the transport equation and the conservation equations for an ar- 
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bitrary choice of the diff usion coefScient i n bot h its dependence on momentum and spatial 



co ordinate was fo und W 



Malkovl fll997h and 



Amato fc Blasil (120051) . 



Amato fc Blasil (120051 ) showed that an excellent approximation to the 



solution f{x,p) has the form 



f{x,p) = fo{p) exp 



Ul 



1-—) I dx 

U2' 



, u{x') 



D{x'^p) 



dhifojp) 



is the local slope of fo{p) in mom entum space. In fact the 



where q{p) -- 

1 — ^ in the argument of the exponential was introduced by 



Blasi. Amato fc Capriolil (120071 ) 



(10) 



actor 



in order to satisfy exactly the boundary condition at the shock, even for weakly modified 
shocks. 

In terms of the distribution function (Eq. [TOj) . we can also write the normalized pressure 
in accelerated particles as: 



An 



Pmax 



dp p^v{p)fo{p) exp 



Pinj 



Xp{x',p) 



(11) 



3D{p,x) 
q(p)uo 



where for simplicity we introduced Xp{x^p) 

By differentiating Eq. [11] with respect to x we obtain 

^ = \{x)^cix)U{x), 
ax 

where 

X{x) =< 1/xp >5, 



I^:r dp P'^)V{p)fo{p) exp [- dx'^^ 



(12) 



(13) 



I^Zr PMp)fo{p) exp [- /° dx'^^^^ 
and U{x) is expressed as a function of C,c{x) through Eq. [21 

Finally, after integration by parts of Eq. [U one is able to express Up{p) in terms of an 
integral involving U{x) alone: 

1 



Upip) 



dx U{x) 



Xp{x,p) 



exp 



"dx'-^ 
Xp{x',p) 



(14) 



which allows one to easily calculate fo{p) through Eq. [3l 

Eqs. [2] and [12] can be solved by iteration in the following way: for a fixed value of the 
compression factor at the subshock, Rsub, the value of the dimensionless velocity at the shock 
is calculated as f/(0) = Rgub/Rtot- The corresponding pressure in the form of accelerated 



1 ( Raub \ 
^ V Rtot ) 



-"la 



. This is used as a 



particles is given by Eq. [2] as ^c(O) = 1 + - -j^^^ 
boundary condition for Eq. [121 where the functions U{x) and A(x) (and therefore fo{p)) on 
the right hand side at the k*^^ step of iteration are taken as the functions at the step (A; — 1). 
In this way the solution of Eq. [2] at the step k is simply 
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(15) 



with the correct hmits when x — > and x — oo. At each step of iteration the functions 
U{x), fo{p), A(x) are recalculated (through Eq. [2], Eqs. [T^ and [3l and Eq. [131 respectively), 
until convergence is reached. The solution of this set of equations, however, is also a solution 
of our physical problem only if the pressure in the form of accelerated particles as given 
by Eq. [2] coincides with that calculated by using the final fo{p) in Eq. [6l This occurs for a 
value of Rsub, which fully determines the solution of our problem for an arbitrary diffusion 
coefficient as a function of location and momentum. 



2.1 Spectra and velocity profiles 

All approaches to particle acceleration at modified shocks predict the formation of a pre- 
cursor in the upstream region, resulting in a gradient of the velocity profile of the fiuid. 
Since qualitatively the spectrum of the accelerated particles is still determined by an effec- 
tive compression factor felt by the particles of given momentum, and the velocity in the 
precursor increases with the distance from the shock, it is easy to infer that the spectrum 
of the accelerated particles is not expected to be a power law and more precisely that it 
should be concave (steeper at low energies and flatter at high energies). Here we discuss the 
detailed shape of the spectrum at the shock as obtained through the two kinetic approaches 
described above. In Fig. [1] (left panel) we plot the spectra as a function of the momentum of 
particles for model B (solid lines) and for model A (dashed lines). The curves are obtained 
for Pmax = lO^rripC, uq = 5 x lO^cm s~^, ^ = 3.5 and for the values of the Mach number 
at upstream inflnity Mq = 10, 100, 1000 (curves labeled as 1, 2 and 3, respectively). Bohm 
diffusion is assumed. The agreement between the two sets of curves is excellent for relatively 
low Mach numbers (Mq ~ 10) and remains good even up to much larger Mach numbers, and 
in fact for all values we have tried. The largest discrepancies between the two methods are 
at the level of ~ 20%. The reason for such discrepancies is to be found in the assumption 
that Upipmax) is required to equal unity in the approach of Model A. 

The velocity proflle of the fluid in the precursor is plotted in Fig. [1] (right panel) for the 
two models (again solid and dashed lines respectively). 

On the X-axis we plot the distance x from the shock in the upstream region in units of 
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10-2 IQO 102 1Q4 lQ-6 10-^ 10-2 102 

p/mc x/xp. 



Figure 1. Left Panel: Spectra of accelerated particles for pmax = lO^mpC, no = 5 X lO^cm s~^, ^ = 3.5 and for Mq = 
10, 100, 1000 (curves l abeled as 1, 2 and 3 r espectively). Bohm diffusion coefficient is adopted. The solid lines are obtained 
with the calculation of I Amato fc BlasH |200^ (model B), the dashed ones with that of Blasi (2002) (model A). Right Panel- 
Velocity profiles in the precursor for the cases in the left panel. Again the solid curves are computed with Model B and the 
dashed ones with Model A. The spatial coordinate is in units of Xpt, with \xpt,\ = D{pmax) / {uoUp{pmax)) as discussed in the 
text. 



Xp*, which is defined as 

I Di^Pmax) 

UoUp{p 

max ) 

Some comments are required on the calculation of U{x) for Model A. As discussed in the 
previous section, this model does not keep any information about the spatial dependence 
of the quantities in the precursor, although such information is somehow contained in the 
relation between a momentum p and the mean diffusion length of particles with such mo- 
mentum, \x{p)\ ~ D{p)/up{p). The dashed lines in Fig. [1] (right panel) are obtained in the 
following way: for a given location x upstream, the equation x = D{p)/up{p) is inverted 
and a corresponding value p of the minimum momentum of particles that may have diffused 
to the point x is obtained. At this point the velocity U{x) (in units of Uq) is by definition 
Up{p) for the value of p corresponding to x. By definition the fiuid velocity in the simple 
model is bound to be unity at x/xp = 1 because no particles are supposed to be able to 
reach farther regions. In the exact solution there is a spread in the distances that can be 
diffusively reached at given momentum and the transition to U{x) = 1 is smoother. This 
difference in the velocity profile affects mainly the results for the spectrum at p ~ Pmax, but 
since these particles carry an appreciable amount of energy in the case of modified shocks, 
the whole spectral shape is somewhat affected (at the level of at most ~ 20% in the strongly 
modified cases). 

As discussed above, the ingredient of model A that makes it much faster in terms of com- 
putational time consists in the assumption that particles with given momentum all travel to 
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some given maximum di stance from th e shock in the upstream region. This physical ansatz 
was first put forward by lEichlerl (jl979l ). This assumption has the precious imphcation that 
the calculation of the spectrum of accelerated particles, of the velocity profile in the precur- 
sor, and of all thermodynamical properties of the background plasma, can be carried out 
without knowing a priori the diffusion coefficient as a function of momentum, which in gen- 
eral is an input to the problem and is very poorly known. The solid curves plotted in Figs. [1] 
are obtained for a Bohm diffusion coefficient in Model B, and as stressed above the agree- 



Blasi. Gabici &: Vannoni 



ment between the two approaches is very good. As discussed by 
( I2OO5I ) the goodness of the result becomes gradually worse for diffusion coefficients in the 
form of Kraichnan {D(p) oc p^/^) and Kolmogorov {D(p) oc p^^^) that show a weaker depen- 



dence on momentum. On the other ha nd, X-ray o 



3servations of SNRs shocks seem to hint 



to Bohm-like diffusion coefficients (see 



Stage et al 



(I2OO6I ) for recent results). 



The price to pay for the short computational time is that one can keep track of the spa- 
tial dependence of the different quantities involved only through the approximate relation 
X ~ —D{p)/up{p). In Model B, which is more challenging computationally, the spatial de- 
pendence is fully accounted for and any form of the diffusion coefficient can be adopted. 
This appears t o be e specially important since the diffusion coefficient as calculated by 



Amato fc Blasil ( 2006 ) in the case of classical streaming instability leading to strong mag- 



netic field amplification has a flat dependence on momentum in the energy region close to 
the maximum momentum (although the overall shape and normalization of the diffusion 
coefficient are close to Bohm-like). 



2.2 Compression factors and multiple solutions 

The stratification of the cosmic ray pressure in the upstream region, together with the 
escape of particles with momentum Pmax from upstream infinity (equivalent to having a 
radiative shock), make the fluid more compressible and therefore lead to an increase of the 
total compression factor Rfot between upstream infinity and downstream. In this section 
we investigate the behavior of Rtot as a function of some parameters of the problem, most 
notably the parameter ^ which defines the injection momentum and the Mach number Mq. 

The fraction of particles that are injected and take part in the acceleration process 
increases when ^ decreases and as a consequence the compression factor also increases. In 
Fig. [2] (left panel) we plot the total compression factor as obtained in Model A. The same 
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Figure 2. Left Panel: Total compression factor as a function of the parameter ^ as obtained with the approach o 
The curves are labelled by the value of log-iQ{pmax /'mpc) adopted. Where multiple solutions are present, the solid line traces 
the most modified one, the dotted line refers to the one close to linear, and t he dashed line to the i ntermediate one (see text 
for discussion). Right Panel: The same curves as obtained with the approach of lAmato Sz BlasH ||2005| ). Both plots are obtained 
for Mo = 100, lio = 5 X 10* cm s"^ 



quantity for the Model B is plotted in the right panel of Fig. O Both plots are obtained for 
Mq = 100, Mo = 5 X lO^cm s~^ and Pmax between 10^ m^c and 10'' rripC (different curves are 
labeled by the value of logiQ^Pmax/iTT'pC))- 

The most striking feature in the behaviour of Rtot versus is the rather sharp transition 
from strongly modified shocks (low values of ^, Rtot ^ 4) to weakly modified shocks (large 
values of ^, Rtot ~ 4). The even more striking point is that, for Pmax > lO^mpC, in the very 
thin transition region in parameter spa ce three soluti o ns rnay app ear. Multiple solutions 
were initi ally found in two-fluid models iDrurv fc Volkl (119801 . Il98ll ) and later found in the 



models of 



Malkovl fll997h : iMalkov. Diamond fc Volkl (j2000|). After introducing the thermal 



leakage model for injection the appearance of multiple solutions was found to be strikingly 



reduced and limited to very narrow regions in parameter space (IBlasi. Gabici &: Vannoni 



( 120051 )). The multiple solutions are here shown to exist even in Model B, though, again, in 
very narrow regions of parameter space. In fact this region is so small that the phenomenon 
was previously missed, since for standard values of the parameters (in particular for ^ = 3.5 
which is most often used) we never find more than one solution. 

Although we cannot prove it in a formal way at the present stage, it is likely that the ap- 
pearance of the multiple solutions in very narrow regions of parameter space is accompanied 
by the existence of some type of instability that allows the system to chose among the three: 
from Fig. [2] one can see that the three solutions are found only in the region of values of ^ 
for which the behaviour of the system suffers a transition from a strongly modified shock 
< C,*) to a weakly modified shock (^ > ^*). In the transition region, a tiny change in the 
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Figure 3. Total compression factor fo r the case at cons tant temperature (Tq = 10'' K) as a function of the Mach number. The 
left and right panel show the results of iBlasil {iooj and lAmato fc Blasil ||2005|) respectively. The two cases J = 3.5 and ^ = 3.8 
are shown. For § = 3.8 multiple solutions appear. In this case the solid curve refers to the most modified one, the dotted curve 
to the one closer to linear and the dashed curve to the intermediate one. 

value of ^ around ^* leads to either one or the other regime. It is therefore likely that at least 
one of the solutions in this transition region may be unstable. Of the three solutions that 
we find, one is certainly non physical (the intermediate solution in Fig. [2]) since it predicts 
that Rtot increases with increasing ^. 

The question of the stability of the other two solutions is sti ll open . Investigat i ons o n 



Toptygin 



fll999h : 



the stability of thes e shoc ks were carried out by iMond fc Druryl (119981 ) ; 
Kang. Jones fc RvJ hoom . 

The behaviour of Rtot as a function of the Mach number is more subtle in that one can 
change the Mach number by fixing the temperature and varying the shock velocity or by 
fixing the shock velocity and changing the temperature of the background gas, although in 
terms of astrophysical applications the first case is probably the most relevant or at least 
the most frequent. Below we consider the two situations separately. 

Increasing the Mach number (for a fixed temperature of the background gas) leads to 
an increasing modification of the shock and therefore to an increase in the value of Rtot- 



For To = 10^ K, p„ 



lO^mpC and ^ = 3.8, the parameters used in the lower curves of 



Fig. [3l this trend continues up to Mq ~ 500. For larger values of the Mach number, there 
is no energy left to convert into accelerated particles and the shock returns to be a test 
particle accelerator. In the thin transition region between the strongly modified regime and 
the regime of weakly modified shocks again three solutions appear. If = 3.5 is used instead 
of = 3.8 the transition moves to much larger Mach numbers, of no astrophysical interest. 
These characteristics are shown in Fig. [3] for the model A (left panel) and B (right panel). 
The two sets of curves are in very good agreement, although for the same parameters model 
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Figure 4. Total compression factor for the c ase at constant v elocity (up = 5 X 10^ cm s~^) as a function of tlic Mach number. 
The left and right panel show the results of iBlasil ll2002h and lAmato fc Blasil ||2005| ) respectively. The two cases f = 3.5 and 
5 = 3.8 are shown. For ^ = 3.8 multiple solutions appear. In this case the solid curve refers to the most modified one, the 
dotted curve to the one closer to linear and the dashed curve to the intermediate one. 



B shows the appearance of the multiple solutions at slightly larger Mach numbers than the 
simple m odel A. The maximurn value of the compression factor is also slightly larger in the 



model of 



Amato &: Blasil (120051 ) than it is in the simple model. The upper curve in both 



plots refers to ^ = 3.5 and shows that no multiple solutions are found in this case. 

When the increase of the Mach number is achieved by fixing the shock velocity and 
changing the temperature, the behaviour of the total compression factor as a function of 
Mo is as shown in Fig. |H For ^ = 3.8, oddly enough, the multiple solutions appear for large 
values of Mq and remain three irrespective of how large Mq becomes. At Mq ~ 100 — 200 
(for the values of the parameters adopted here) there is a bifurcation: one branch that 
smoothly connects to the weakly modified solution for low values of Mq remains and tends 
asymptotically to Rtot ~ 5, while two other branches appear, one with compression factor 
that keeps increasing and the other with compression factor that tends to Rtot ~ 15- 

Even a qualitative comparison of Figs. [3] and H] reveals that the two cases are intrinsically 
different. For C, = 3.5 (upper curve) no multiple solutions are found, and the shock always 
shows an appreciable level of modification. 

The presence or absence of a transition to a weakly modified shock at sufficiently large 
Mach number can be understood in a semi- quantitative way by using the following argument. 
Let us assume that the shock is weakly modified, namely that the pressure in the form of 
accelerated particles is much smaller than PqUq, so that the spectrum can be approximated 
as a power law with slope ~ 4 (/o(p) ~ P~'^)- In this case the total pressure of the accelerated 
particles at the shock location is 



CR 



fPmax 



Pinj 



dpp^v(p)fo{p) 
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Pmax 



(16) 



The term 47rp?„j/o(pj„j) is roughly the total number of particles at the shock location, which 
can be written as ~ r]niUi/u2 ~ Arjui. It follows that 



Pgr ~ -c4:r]nipinj In 



Pr, 



nripC ^ 



(17) 



If one neglects the logarithmic term, the pressure in accelerated particles is easily seen to 
scale with the injection momentum only. In the context of the thermal leakage approach 
and in the limit of strong shock, we can write p^j = ^j2mpkT2, where kT2 = (3/16)Mp 



It follows that Pinj = ^y(3/8)mpMi, so that Pcr oc Ui. The linear scaling of the cosmic 
ray pressure with ui should be compared with the ram pressure scaling, which is piuf. 
This comparison immediately suggests that for sufficiently large values of Ui there is always 
enough kinetic pressure to fuel cosmic ray acceleration without appreciable modification of 
the shock. Notice that when the Mach number is increased by keeping the fiuid velocity 
constant, this argument does not apply and indeed in that case the transition is not seen, 
while the three solutions persist at sufficiently high Mach numbers. If the Mach number 
is increased by increasing the shock velocity (keeping the temperature constant) then the 
transition is observed when the fluid velocity is larger (in fact about one order of magnitude 
larger) than a critial value u^, which can be estimated from Eq. [TT) 



c 



Pr, 



(18) 



It is clear that the larger the values of ^ the smaller the critical velocity for which the 



shock may become weak. 



Berezhko k Ellison! (119991 ). 



y modifled. The same threshold effect was previously found by 



Some more discussion is required about the role of stationarity in the appearance of 



multiple solutions. The p henomenon of multip 



of both two fluid models 



Drurv fc Volkl (11980 



e solutions has been reported in the context 



198ll ) and kinetic approaches, but in both 



cases stationarity was assumed. Multiple solutions do not seem to appear in time- dependent 
numerical calculations. From the physical point of view a stationary solution of the diffusion- 
convection equation when particle escape or energy losses are absent cannot exist. Even in 
the context of the test-particle approach quasi-stationarity can be recovered only at low 
momenta, far from Pmax-, while the temporal evolution reflects into an increase of Pmax- In the 
non-linear regime an increase of Pmax leads to a modiflcation of the precursor and therefore of 
the entire spectrum. Hence the assumption of stationarity is harder to justify in the context 
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of a non-linear theory of particle acceleration at shocks. This internal inconsistency of both 
two-fluid models and semi-analytical kinetic models may be one of the reasons why multiple 
solutions do not appear in time-dependent approaches. 



2.3 Advected and escaping fluxes 

One of the most important predictions of the non-linear theory of diffusive particle accelera- 
tion at shocks is that for strongly modified shocks an appreciable fraction of energy is in the 
form of particles at the maximum momentum. On the other hand, these particles are also 
the only ones that are allowed to leave the system from upstream infinity. In other words, 
the return probability to the shock from upstream is unity for all particles but for those with 
P ~ Pmax- The energy carried away from the shock by the highest energy particles makes 
the shock radiative and the increased compressibility of the background gas enhances the 
modification of the shock structure. The equation for the conservation of the flux of energy 
provides us with precious information, namely the flux of energy in the form of accelerated 
particles that is advected towards downstream infinity and the one that escapes towards 
upstream infinity. 

The conservation equation between downstream and upstream infinity can be written in 
the following form: 

^p2ul H '^^Pg,2U2 H '^^Pc,2U2 = ^PO^iQ -^^9,0^0 " Fe, (19) 

2 7s - 1 7c - 1 2 73-1 

where Fe is the flux of particles escaping at the maximum momentum from the upstream 
section of the fluid (Berezhko & Ellison, 1999). This term is peculiar of modified shocks, 
being completely negligible when acceleration takes place in the test particle regime, as we 
confirm in the calculations below. 

In Eq. [in] we can divide all terms by (l/2)po^o ^^"^ calculate the normalized escaping 
flux: 

^' = 2 _ J 1^£eA _ 1 l^E^ (20) 

^ Rm M^ilg - 1) Rtot Ig - 1 PoMo Rtot 7c " 1 Po^o ' 

From momentum conservation at the subshock we also have: 

Pc,2 _ Rsub _ 1 _|_ 1 f Rsub\ '^^ 

Poul Rtot Rtot IgM^ V Rtot ) 

SO that the escaping flux only depends upon the environment parameters (for instance the 
Mach number at upstream infinity) and the compression parameter Rsub which is part of 
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the solution. The adiabatic index appropriate for cosmic rays, 7c, is here calculated self- 
consistently as: 

^ . , ^ . , lg:rdpA7rpMp)foip) 

where Ec is the energy density in the form of accelerated particles and e{p) is the kinetic 
energy of a particle with momentum p. It can be easily seen that 7^ -—>■ 4/3 when the energy 
budget is dominated by the particles with p ~ Pmax (namely for strongly modified shocks) 
and 7c ^ 5/3 for weakly modified shocks. In Eq. [20] the term F^^^ = -^-^y^^ou^ clearly 
the fraction of flux which is advected downstream with the fluid. 

The escaping flux, the advected flux and the total flux in the form of accelerated particles 
are plotted in Fig. O for mq = 5 x lO^cm s~^, Pmax = 10® mpC and ^ = 3.5. All fluxes are 
normalized to {l/2)poul. Clearly the difference between the total flux in accelerated particles 
and unity gives the rate of conversion of the total energy into heating of the background 
gas. For large Mach numbers this difference vanishes, which implies that the gas is not 
appreciably heated at the shock, one of the most impressive predictions of the non-linear 
theory of particle acceleration. 

It is worth discussing in some details the physical meaning of the advected and escaping 
fluxes. As we already pointed out earlier, the escaping flux is all in the form of a narrow 
distribution in momentum around p ~ Pmax- One should keep in mind that in situations 
of astrophysical interest, such as in the case of Supernova Remnants (SNRs), after the 
beginning of the Sedov phase the maximum momentum decreases with time. A distant 
observer is likely to observe (from a single SNR) an overlap of p eaked functions with p ~ 
Pmax{t)- In specific situations (e.g. iPtuskin &: Zirakashvilil (120051 )) it has been shown that 
this superposition leads to power law time-integrated spectra. Since the particles of p ~ Pmax 
carry an appreciable fraction of energy (in Fig. O one can see that —>■ 1 for large Mach 
numbers), this spectrum should roughly coincide with the injection spectrum emitted by the 
SNR, despite the fact that the spectrum at the shock has a shape of the type illustrated in 
Fig. [1] (at a given time). One should keep this argument in mind when arguing about the 
effects induced by concavity in the source spectra on the spectrum of diffuse cosmic rays 
observed at the Earth. 

The spectrum of particles advected downstream, roughly speaking is ~ U2fo{p), therefore 
the flux of energy carried by these particles is ~ U2-Pc = -^Pc^o- Since for modified shocks 
Rtot ^ 1, the actual advected flux is not very large. We stress that this is the flux which is 
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actually useful from the phenomenological point of view in terms of conversion of energy into 
gamma rays and other types of radiation within the source. On the other hand the particles 
that escape from upstream may i nteract with a thick target i n the upstream region an d 

mom . 



Moskalenko et al. 



produce a detectable signal there (IGabici fc AharonianI (120071 ): 
From Fig.[5]one can see that the advected flux takes less than ~ 10% of the total energy influx 
for strongly modified shocks. Actually such flux is larger for shocks which are less modified, 
which in Fig. [5] correspond to lower Mach numbers, but even in this case it remains smaller 
than ~ 30 - 40% in units of (l/2)po^^o. 



The particles which are advected downstream remain behind the shock and may even- 
tually leave the remnant only at later times, suffering adiabatic energy losses due to the 
expansion of the supernova shell. Their final spectrum, as could be observed by a distant 
observer, is a complex convolution of the temporal evolution of the shell, the changing max- 
imum momentum, and the different levels of shock modification at different times. 



For low values of the Mach number the shock modification predicted by kinetic models 
decreases, and as a consequence the concavity typical of the modified spectra is reduced, 
until the test particle solution is approached. Once the slope of the spectrum at p ~ Pmax 
drops below ~ 4 the energy flux that escapes towards upstream infinity gets suppressed, as 
shown in Fig. [H 



The physical origin of the escaping flux is again rather puzzling and requires some com- 
ments: from the mathematical point of view the requirement of having an escaping flux at 
upstream infinity derives from imposing energy conservation. From the physical point of 
view it is not obvious that the conditions for such escape exist. For instance, in the case of 
supernova explosion, one may envision escape of particles during the Sedov phase, but not 
during the free expansion phase. In fact in the latter phase, the maximum momentum in- 
creases with time. Nevertheless, a time-independent approach such as the stationary kinetic 
models discussed here, would predict a flux escape even during the free expansion phase, in 
order to conserve energy. Although this problem is usually not discussed in the literature, 
we think that this may again call for the need to develop time-dependent non-linear kinetic 
models. In such approaches particle escape from upstream infinity should be invoked only 
when there are the physical conditions for it to occur. 
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Figure 5. Normalized total (solid line), advected (dashed line) and escaping (dash-dotted line) flux of accelerated particles as 
functions of the Mach number. 



2.4 The effect of turbulent heating on the energy fluxes 

In the previous section we have amply discussed that the non-linear theory of particle acceler- 
ation tends to convert large fractions of the flux traversing a shock into accelerated particles. 
One of the effects that in astrophysical situations are likely to reduce such an efficiency is 
the so-called turbulent heating. This generic expression is used to refer to any process that 
may determine non-adiabatic gas heating in the precursor. The tw o best known examples of 
this type of pro c esses are Alfven heating (IMcKenzie &: Volkl (119821 )) and acoustic instability 
( iDrury fc Falld (119861 )). Both effects are however very hard to implement in a quantitative 
calculation: in the case of Alfven heating, the mechanism was originally introduced as a 
way to avoid the turbulent magnetic fleld to grow to non-linear levels, while it is usually 
used even in those cases in which 6B/Bq ^ 1. Acoustic instability develops in the pressure 
gradient induced by cosmic rays in the precurs or and results in the development of a train 
of shock waves that heat the background gas (IDrury Sz Falld (119861 )). The analysis of the 
instability is carried out in the linear regime, therefore it is not easy to describe quantita- 
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tively the heating effect. In both cases the net effect is the non-adiabatic heating of the gas 
in the precursor, which resuhs in the weakening of the precursor itself and in the reduction 
of the acceleration effciency compared with the case in which the turbulent heating is not 
taken into account. 

In order to simply illustrate the effect and compare the findings of model A and B in 
describing the fluxes of advected and escaping a ccelerated particles i n the case of turbulent 



heating, we adopt the simple recipe provided in iBerezhko fc EUisonl (119991 ) for Alfven heat- 
ing. The recipe consists in modifying the relation between the gas pressure upstream of the 
subshock at the location x, Pg{x), and the gas pressure at upstream infinity, in order to take 



into account the amount of non adiabatic heating. The proposed expression is 



9,0 \ Po J [ ^vIa,0 



7s' 



1 - 



(23) 



.Pi^)/ 

where Ma,o is the Alfvenic Mach number at upstream infinity. One can easily check that the 
non-adiabatic heating vanishes for Ma,o — > oo. With this simple recipe, Eq. [2] is modified in 
the following way: 

If one keeps in mind that this modification also changes the relation between and Rtot, 
it is easy to realize that the computational procedures of the two kinetic models are left 
otherwise unchanged and we are now able to determine the effect of turbulent heating, at 
least in the context of this simple approach. 

The best way to illustrate the effect of turbulent heating is by plotting the same fluxes 
as in Fig. [5] but including now the turbulent heating. The results are plotted in Fig. [6|, where 
the lines are obtained with Model A and the symbols represent the preditions of Model B. 
Aside from the very good agreement between the two methods, it is worth stressing that 
the turbulent heating leads to less modified shocks, and therefore to a smaller flux towards 
upstream infinity. The total flux (advected plus escaping) in the form of accelerated particles 
is also substantially reduced, so that appreciable heating can occur at the subshock. 



3 THE MAXIMUM MOMENTUM OF THE ACCELERATED PARTICLES 



The maximum mo mentum that particles accelera ted at modified shocks can achieve has 



been calculated by 



Blasi. Amato fc Capriolil (120071 ). for arbitrary level of modification and 



for arbitrary choice of the diffusion coefficient. Here we summarize the calculations of the 
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Figure 6. Normalized total (solid line), advected (dashed line) and escaping (dash-dotted line) flux of accelerated particles 
as functions of the Mach number when turbulent heating is tak e n into account. The assumed value of the magnetic field is 
Bp = 10/jG. The line s are obtained with the approach of iBlasil (|2002|) . while the symbols are the results of the method of 
lAmato fc Blasil ||2005| '). 



acceleration time and illustrate a recipe that allows us to include this calculation in Model 
A. We compare the results to those obtained with the formal approach of Model B. The 
acceleration time up to a momentum p is given by 



< t >= - 



dh 



ds 



3R 



tot 



RtotD2{p') 



+ 



uoA{p') 



fp dp' j 

^=0 "-u ■Jp^nJ P' \ RtotU.p{p') - 1 ' RtotUpip') - 1 ^ 
where D2{p) is the diffusion coefficient in the downstream plasma and 



dxexp I (l 
[ 3 



U2 
Ui 



(25) 



(26) 



10 D{x') J 

In the limit of test particle acceleratio n, this accelerat ion time reduces to the well known 
expression ( Lagage fc Cesarsky ( IQSsJ lbl): Iprury ( 1983 )): 



< t >-- 



U2 Ui . 



Ui - U2 

We recall that in the model of 



(27) 



Blasil (120021 ) the information on the spatial distributions in 
the precuror is kept only through the relation x{p) = D{p) / (uqUp). This recipe has also been 
used to couple x to p in the integral for A{p). 
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The maximum momentum is assumed here to be determined by the equahty between the 
acceleration time and the age of the accelerator. Clearly other recipes can be easily imple- 
mented as well. Below we consider two cases: 1) a background magnetic field Bq = lO/iG; 2) 
the fie d is amplified by stream i ng ins tability to the saturation value SB = Bq [2MA_fl-^^ 



flBelll fll978af ): lAmato fc Blasil fl2006h l 



In both cases the diffusion coefficient is assumed to be Bohm-like. The magnetic field 
downstream, needed to calculate D2{p), remains Bq in case 1), while in case 2) the two 
components of the field perpendicular to the shock normal are compressed by Rsub whereas 
the parallel component is left unaltered. 

Our results for the maximum momentum as a function of the Mach number are plotted in 
Fig. [71 where we assumed that the diffusion coefficient is spatially constant in the precursor. 
We also adopted uq = 5x lO^cm s^^, ^ = 3.5 and Bq = 10 fiG, while the age of the accelerator 
is fixed at 10 00 years. The curves in Fig. [7| r efer to Model A, while the symbols illustrate 



the results of 



Blasi. Amato fc Capriolil (120071 ) for Model B. The lower curve (and symbols) 



refers to the case Bq = lO/iG, while the upper curve (and symbols) refers to the case of 
magnetic field amplified by streaming instability. Once again, the agreement between the 
two kinetic approaches is very good, once Model A is completed with a recipe that allows 
us to recover the spatial information in the precursor. 



As already noticed by iBlasi. Amato fc Capriolil (120071 ). the main factor in enhancing 



the maximum momentum of the accelerated particles is the magnetic field amplification. 
On the other hand, increasing the Pmax also induces a larger shock modification (the same 
effect occurs when the Mach number increases). Larger shock modifications lead to a slight 
reduction of the maximum momentum. The reason is simple to understand: as an order of 
magnitude estimate, the acceleration time is inversely proportional to the square of a mean 
fiuid velocity, as averaged over the precursor. For strongly modified shocks the fiuid upstream 
slows down appreciably, thereby increasing the acceleration time. For our benchmark values 
of the parameters, typical of supernova remnants, we obtain maximum momenta between 
5 X 10^ and 2 x 10^ rripC, comparable with the energy where the knee is observed in the 
cosmic ray spectrum. 
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Figure 7. Maxi mum momen tum of the accelerated particles for spat i ally c onstant diffusion coefficient. The lines show the result 
of the method of Blasi| | |2002|) . the symbols refer to llAmato fc Bias] l|2005)). The lower line (and symbols) refer to Bq = lOfiG, 
while the upper curve (and symbols) are obtained for magnetic field amplified by streaming instability. 



4 CONCLUSIONS 

The recent growth of observational evidence of efficient particle acceleration, at least in the 
case of supernova remnants, has posed a serious challenge to describe the chain of physical 
processes that are all taking place at the same time and in a strongly correlated manner: 
particle acceleration is likely to excite streaming instability leading to magnetic field ampli- 
fication upstream. This field is then advected downstream and leads to effective synchrotron 
emission of electrons in narrow filaments which are observed in X-rays. The amplified mag- 
netic field makes it possible to reach maximum momenta of accelerated protons which are 
comparable with the knee in the cosmic ray spectrum, but it also leads to, and requires, 
efficient particle acceleration, which in turn modifies the shock structure. The description of 
this complex chain of non-linear phenomena ideally requires the particle acceleration process 
to be treated at the same time of and coupled to a hydrodynamic (or even MHD) code for 
the evolution of the remnant, thereby leading to the need for a fast, efficient and accurate 
code for the description of non-linear particle acceleration. In many previous papers the 
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simple approach of iBerezhko fc EUisonl (119991 ) was adopted. That approach, though quah- 
tatively appropriate, forced the spectrum of accelerated particles to be a broken power law 
with given points where the slopes changed, an assumption that can only be considered as 
a working hypothesis. 
The approach of 



Blasil ( 2002 



2004 ). discussed here as Model A, is based on a physical 



ansatz on the spatial distribution of particles in the precursor. It allows one to obtain the 
spectrum of accelerated particles and all thermodynamical quantities of the fluid in a short 
comput ational time. For thi s reason it has recently been implemented in the hydrodynamical 



code of 



Ellison et al. 



(120071 ). Model B, also discussed here, leads to a formal solution of the 
problem of particle acceleration, that we have shown to work well for a number of different 
assumptions on the spatial and momentum dependence of the diffusion coefficient. However 
its computational time is too long to allow one to use it in more complex calculations. 

In this paper we discussed the main results obtained by using the two kinetic approaches, 
with two goals in mind: 1) show that Model A provides sufficiently accurate results to allow 
its use in more complex calculations; 2) investigate one aspect of kinetic approaches (and in 
fact of all stationary approaches) that is still poorly understood, namely the appearance of 
multiple solutions. 

In ord er to assess the goodness o f Model A we compared its results with the fomal 



solution of 



Amato fc Blasil ( 



200 



2 



20061 ) (Model B). This was possible after completing Model 



A with a recipe to infer the spatial information on particle distribution in the precursor, 
as discussed in §2.11 The differences between the results of the two models are typically 
smaller than 20% for all quantities which have been calculated (spectra, velocity profile in 
the precursor, compression factors). The recipe mentioned above a lso allowed us to apply to 
Model A the calculation of the acceleration time as presented by iBlasi. Amato fc Caprioli 
(120071 ) and therefore to infer the maximum energy of the accelerated particles. Also in this 
respect, as discussed in ^ Model A returns results in very good agreement with the formal 
solution (Model B). 

The most ineteresting insights came however from the investigations on the presence of 
multiple solutions and escaping fluxes. 



M ultiple solutions were first found in the conte xt of stat i onary two fluid models 



19801 . 



(1200 



Drurv fc Volk 



Malkov. Diamond fc Volk 



■9811) and later in the kinetic ap proach of iMalkovl (119971 ): 
Blasi. Gabici fc Vannoni (2W) showed that treating injection as a thermal leakage. 



multiple solutions persisted only in a very narrow range of parameters. While confirming 
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this finding here for Model A, we found that also the formal solution of Model B leads 
to multiple solutions, in about the same regions of parameters as for Model A. Since the 
two methods are based on quite different criteria for the solution of the equations, this is 
a strong hint to the fact that the multiple solutions are an intrinsic property of the system 
and not an artifact of the iteration procedure used to solve the equations. A detailed discus- 
sion of the appearance of multiple solutions has been presented in §2.21 However a physical 
understanding of the multiple solutions is still missing. It is intriguing that they are found 
in time-independent approaches, where stationarity is assumed, while they are not found 
in time-dependent approaches. One should keep in mind that a stationary solution of the 
equations for particle acceleration at a shock, in the absence of losses and escape does not 
exist, and that the assumption of stationarity is therefore rather artificial. 

On the other hand, time-dependent approaches are usually based on the solution of the 
coupled transport equation and fluid equations in such a way that the solution at a given 
time t is advanced, following a predefined integration scheme, to a time t + At. In such 
approaches we cannot envision a procedure that would lead to the appearance of multiple 
solutions. What probably could happen is that there may be a strong dependence on initial 
conditions. To our knowledge there has been no investigation of such effects. Moreover, as 
pointed out above, even if multiple solutions do exist one or more of them are likely to be 
unstable. 

The requirement of stationarity implies that a flux of energy must escape from upstream 
infinity. Such a flux is actually predicted on physical grounds in some circumstances, such 
as the slowing down of the fluid motion in the Sedov phase of a SNR evolution, or because 
of the fact that during such phase the magnetic field amplification is expected to decrease 
with time, thereby reducing Pmax- In this way particles accelerated to higher maximum 
momentum at previous times can no longer be confined in the accelerator and escape. But 
the stationary non-linear calculations of particle acceleration at modified shocks also predict 
escape in situations where it is not immediate to foresee it on physical grounds. Again, since 
the role of this escaping flux has profound implications on the acceleration process it would 
be appropriate to investigate it using time-dependent techniques. 
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